INLA setup

hyper.iid = list(theta = list(prior = "pc.prec", param = c(1, 0.01)))
hyper.bym2 = list(theta1 = list("PCprior", c(1, 0.01)), 
                  theta2 = list("PCprior", c(0.5, 0.5)))
threads = parallel::detectCores()

Model 1.0: no covariates

We use a model structure similar to Poisson regression, where \(O_{i,j,k}\), the number of observed deaths in 2020 in municipality \(i\), age group \(j\) and sex group \(k\), depends on the number of expected deaths \(E_{i,j,k}\) based on historical data and a linear predictor \(\lambda\).

\[ O_i \sim \text{Poisson}(E_i \times \exp(\lambda)) \\ \] At start, the linear predictor \(\lambda\) only includes one intercept parameter \(\alpha\), so that the estimate of \(\exp(\alpha)\) can be interpreted as an average relative excess mortality for 2020. By adding covariates to \(\lambda\), we aim to disentangle the various factors that are associated with excess mortality at the local level.

We implement this model in R-INLA, a Bayesian inference package that is especially adapted to spatial data. This is achieved in practice by including \(\log (E_{i,j,k})\) as an offset (although an alternative formulation based on the E argument exists). During model development, we compare different model versions based on the WAIC (lower values imply a better fit).

model1.0 = INLA::inla(munici_observed ~ 1 + offset(E),
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)
summary(model1.0)

Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, 
   ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
   scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
   ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
   verbose, ", " lincomb = lincomb, selection = selection, control.compute 
   = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, 
   control.fixed = control.fixed, ", " control.mode = control.mode, 
   control.expert = control.expert, ", " control.hazard = control.hazard, 
   control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " 
   control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
   ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
   num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, 
   working.directory = working.directory, ", " silent = silent, inla.mode 
   = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = 
   .parent.frame)") 
Time used:
    Pre = 0.625, Running = 0.463, Post = 0.083, Total = 1.17 
Fixed effects:
             mean    sd 0.025quant 0.5quant 0.975quant  mode kld
(Intercept) 0.303 0.004      0.296    0.303       0.31 0.303   0

Watanabe-Akaike information criterion (WAIC) ...: 95846.69
Effective number of parameters .................: 6.44

Marginal log-Likelihood:  -48505.00 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.0$summary.fixed)

                mean       sd 0.025quant 0.5quant 0.975quant     mode kld
(Intercept) 1.354092 1.003662   1.344425 1.354092   1.363829 1.354092   1

As a sanity check, we find a relative excess mortality of 35% for 2020, that is coherent with a simple calculation (74,776 observed / 55,221 expected = 1.35). Remember that we excluded the age group 0-40, which explains why this is higher than reported numbers in Switzerland, generally around 10% for 2020. We can also look at the model fit and at the residuals. Obviously the model fit is not good here, as this basic model assumes a unique relative excess mortality for all areas, sexes and age groups.

Model 1.1: age and sex

We hypothesize that excess mortality affected different age and sex groups differently. We thus add the age group, the sex and the interaction of the two as covariates.

model1.1 = INLA::inla(munici_observed ~ - 1 + offset(E) +
                        sex:age_group,
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)
summary(model1.1)

Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, 
   ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
   scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
   ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
   verbose, ", " lincomb = lincomb, selection = selection, control.compute 
   = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, 
   control.fixed = control.fixed, ", " control.mode = control.mode, 
   control.expert = control.expert, ", " control.hazard = control.hazard, 
   control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " 
   control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
   ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
   num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, 
   working.directory = working.directory, ", " silent = silent, inla.mode 
   = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = 
   .parent.frame)") 
Time used:
    Pre = 0.35, Running = 0.634, Post = 0.0681, Total = 1.05 
Fixed effects:
                           mean    sd 0.025quant 0.5quant 0.975quant   mode kld
sexFemale:age_group40-59 -0.003 0.024     -0.050   -0.003      0.045 -0.003   0
sexMale:age_group40-59    0.226 0.018      0.190    0.226      0.262  0.226   0
sexFemale:age_group60-69  0.029 0.020     -0.009    0.029      0.068  0.029   0
sexMale:age_group60-69    0.387 0.015      0.358    0.387      0.416  0.387   0
sexFemale:age_group70-79  0.251 0.013      0.226    0.251      0.276  0.251   0
sexMale:age_group70-79    0.392 0.011      0.371    0.392      0.413  0.392   0
sexFemale:age_group80+    0.480 0.006      0.468    0.480      0.492  0.480   0
sexMale:age_group80+      0.140 0.007      0.126    0.140      0.154  0.140   0

Watanabe-Akaike information criterion (WAIC) ...: 94024.83
Effective number of parameters .................: 17.32

Marginal log-Likelihood:  -47606.94 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.1$summary.fixed)

                              mean       sd 0.025quant  0.5quant 0.975quant
sexFemale:age_group40-59 0.9972893 1.024455  0.9511633 0.9972893   1.045652
sexMale:age_group40-59   1.2533073 1.018531  1.2090047 1.2533073   1.299233
sexFemale:age_group60-69 1.0297260 1.019763  0.9909767 1.0297260   1.069990
sexMale:age_group60-69   1.4728023 1.015056  1.4302909 1.4728023   1.516577
sexFemale:age_group70-79 1.2857289 1.012778  1.2541277 1.2857289   1.318126
sexMale:age_group70-79   1.4799349 1.010613  1.4496266 1.4799349   1.510877
sexFemale:age_group80+   1.6164112 1.006044  1.5974333 1.6164112   1.635615
sexMale:age_group80+     1.1503737 1.007045  1.1346542 1.1503737   1.166311
                              mode kld
sexFemale:age_group40-59 0.9972893   1
sexMale:age_group40-59   1.2533073   1
sexFemale:age_group60-69 1.0297260   1
sexMale:age_group60-69   1.4728023   1
sexFemale:age_group70-79 1.2857289   1
sexMale:age_group70-79   1.4799349   1
sexFemale:age_group80+   1.6164112   1
sexMale:age_group80+     1.1503737   1

model1.1$waic$waic - model1.0$waic$waic 

[1] -1821.858

As expected, the relative excess mortality varies a lot across age and sex groups. It’s very small in females aged 40-59 and 60-69 (in fact the data is compatible with no excess in both cases). It increases in females aged 70-79, and even more so aged 80+. It’s comparatively higher than females in males below 80, but somewhat surprisingly lower in age group 80+. We observe an improvement of the model fit, not easy to see on the plot because of the large number of points, but made clear by the large decrease in WAIC.

Model 1.2: spatial variability

We now account for spatial variability, first in a simple way using an i.i.d. random effect, so that all municipalities can vary independently from each other around a global average.

model1.2 = INLA::inla(munici_observed ~ - 1 + offset(E) +
                        sex:age_group +
                        f(id_space, model = "iid", hyper = hyper.iid, constr = TRUE),
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)
summary(model1.2)

Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, 
   ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
   scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
   ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
   verbose, ", " lincomb = lincomb, selection = selection, control.compute 
   = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, 
   control.fixed = control.fixed, ", " control.mode = control.mode, 
   control.expert = control.expert, ", " control.hazard = control.hazard, 
   control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " 
   control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
   ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
   num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, 
   working.directory = working.directory, ", " silent = silent, inla.mode 
   = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = 
   .parent.frame)") 
Time used:
    Pre = 0.398, Running = 1.16, Post = 0.143, Total = 1.7 
Fixed effects:
                          mean    sd 0.025quant 0.5quant 0.975quant  mode kld
sexFemale:age_group40-59 0.001 0.025     -0.048    0.001      0.049 0.001   0
sexMale:age_group40-59   0.234 0.019      0.197    0.234      0.271 0.234   0
sexFemale:age_group60-69 0.034 0.020     -0.005    0.034      0.073 0.034   0
sexMale:age_group60-69   0.398 0.015      0.368    0.398      0.429 0.398   0
sexFemale:age_group70-79 0.256 0.013      0.230    0.256      0.283 0.256   0
sexMale:age_group70-79   0.401 0.011      0.380    0.401      0.423 0.401   0
sexFemale:age_group80+   0.494 0.007      0.480    0.494      0.508 0.494   0
sexMale:age_group80+     0.151 0.008      0.136    0.151      0.167 0.151   0

Random effects:
  Name    Model
    id_space IID model

Model hyperparameters:
                        mean   sd 0.025quant 0.5quant 0.975quant  mode
Precision for id_space 50.33 5.42      40.80    49.95      61.99 49.21

Watanabe-Akaike information criterion (WAIC) ...: 93744.68
Effective number of parameters .................: 812.83

Marginal log-Likelihood:  -47400.34 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.2$summary.fixed)

                             mean       sd 0.025quant 0.5quant 0.975quant
sexFemale:age_group40-59 1.000593 1.024880  0.9535134 1.000593   1.049999
sexMale:age_group40-59   1.264004 1.019011  1.2181738 1.264003   1.311563
sexFemale:age_group60-69 1.034804 1.020225  0.9949571 1.034804   1.076249
sexMale:age_group60-69   1.489276 1.015573  1.4448262 1.489274   1.535106
sexFemale:age_group70-79 1.292334 1.013406  1.2590236 1.292333   1.326532
sexMale:age_group70-79   1.494002 1.011222  1.4616657 1.493999   1.527066
sexFemale:age_group80+   1.638953 1.007199  1.6160785 1.638945   1.662193
sexMale:age_group80+     1.163545 1.007955  1.1456118 1.163542   1.181775
                             mode kld
sexFemale:age_group40-59 1.000593   1
sexMale:age_group40-59   1.264002   1
sexFemale:age_group60-69 1.034804   1
sexMale:age_group60-69   1.489270   1
sexFemale:age_group70-79 1.292332   1
sexMale:age_group70-79   1.493995   1
sexFemale:age_group80+   1.638930   1
sexMale:age_group80+     1.163536   1

model1.2$waic$waic - model1.1$waic$waic 

[1] -280.1541

The age and sex effect remains similar, but the model fit as measured by the WAIC is improved now that we account for local differences. We can observe this municipality effect, that applies in all age and sex groups of the municipality in the same way.

We find noisy estimates in some places, suggesting issues related to small area estimation. One solution is to partially pool information between municipalities that are geographically linked.

Model 1.3: structured spatial variability

We still focus on spatial variability, but now the municipalities are no longer independent: we account for the correlation between neighboring municipalities with a BYM model. This will allow us to differentiate between what can be attributed to a municipality in particular, and what can be attributed to regional effects (like a COVID wave).

model1.3 = INLA::inla(munici_observed ~ - 1 + offset(E) +
                        sex:age_group +
                        f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                          hyper = hyper.bym2, constr=TRUE),
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)
summary(model1.3)

Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, 
   ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
   scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
   ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
   verbose, ", " lincomb = lincomb, selection = selection, control.compute 
   = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, 
   control.fixed = control.fixed, ", " control.mode = control.mode, 
   control.expert = control.expert, ", " control.hazard = control.hazard, 
   control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " 
   control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
   ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
   num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, 
   working.directory = working.directory, ", " silent = silent, inla.mode 
   = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = 
   .parent.frame)") 
Time used:
    Pre = 9.71, Running = 3.29, Post = 0.432, Total = 13.4 
Fixed effects:
                          mean    sd 0.025quant 0.5quant 0.975quant  mode kld
sexFemale:age_group40-59 0.003 0.025     -0.045    0.003      0.051 0.003   0
sexMale:age_group40-59   0.234 0.019      0.197    0.234      0.271 0.234   0
sexFemale:age_group60-69 0.036 0.020     -0.003    0.036      0.076 0.036   0
sexMale:age_group60-69   0.401 0.016      0.370    0.401      0.431 0.401   0
sexFemale:age_group70-79 0.258 0.013      0.232    0.258      0.285 0.258   0
sexMale:age_group70-79   0.403 0.011      0.381    0.403      0.425 0.403   0
sexFemale:age_group80+   0.495 0.007      0.481    0.495      0.510 0.495   0
sexMale:age_group80+     0.155 0.008      0.139    0.155      0.171 0.155   0

Random effects:
  Name    Model
    id_space BYM2 model

Model hyperparameters:
                         mean    sd 0.025quant 0.5quant 0.975quant   mode
Precision for id_space 60.710 8.060     46.635   60.084      78.33 58.718
Phi for id_space        0.549 0.098      0.358    0.549       0.74  0.543

Watanabe-Akaike information criterion (WAIC) ...: 93570.08
Effective number of parameters .................: 621.16

Marginal log-Likelihood:  -46477.81 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.3$summary.fixed)

                             mean       sd 0.025quant 0.5quant 0.975quant
sexFemale:age_group40-59 1.003025 1.024914  0.9557685 1.003025   1.052622
sexMale:age_group40-59   1.263648 1.019078  1.2176755 1.263646   1.311367
sexFemale:age_group60-69 1.036956 1.020277  0.9969271 1.036955   1.078596
sexMale:age_group60-69   1.492998 1.015652  1.4482184 1.492994   1.539184
sexFemale:age_group70-79 1.294762 1.013505  1.2611503 1.294760   1.329284
sexMale:age_group70-79   1.496164 1.011349  1.4634298 1.496159   1.529660
sexFemale:age_group80+   1.640505 1.007407  1.6169780 1.640489   1.664471
sexMale:age_group80+     1.167878 1.008154  1.1494397 1.167872   1.186645
                             mode kld
sexFemale:age_group40-59 1.003024   1
sexMale:age_group40-59   1.263643   1
sexFemale:age_group60-69 1.036953   1
sexMale:age_group60-69   1.492986   1
sexFemale:age_group70-79 1.294755   1
sexMale:age_group70-79   1.496149   1
sexFemale:age_group80+   1.640455   1
sexMale:age_group80+     1.167861   1

model1.3$waic$waic - model1.2$waic$waic 

[1] -174.6029

We see that the structure accounts for about half of the spatial variability (Phi estimated to 0.5). This addition also improves the model fit as measured by the WAIC.

We observe that many of the municipalities with higher relative excess mortality are in the western and southern parts, the ones that were hit first by COVID-19 in 2020. We also observe areas with higher excess in the North and Northeastern parts. These largely correspond to areas that were hit the most during the first and the second COVID-19 waves of spring and fall 2020 (REF Konstantinoudis et al NatCom 2022).

Model 1.4: local characteristics

Having accounted for regional variability (arguably caused by COVID-19 waves of different timings and scales), we move on to explore the effect of local characteristics at the municipality level. We look at the following covariates (at the municipal level):

• population density in three classes (high, medium or low population density, reference is medium) as defined for each municipality by the Federal Statistical Office;

• median socio-economic position (highest means higher position) in quintiles (reference is 3rd quintile);

• municipality bordering another country.

model1.4 = INLA::inla(munici_observed ~ - 1 + offset(E) +
                        sex:age_group +
                        f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                          hyper = hyper.bym2, constr=TRUE) +
                        density_high + density_low +
                        sep1+sep2+sep4+sep5 +
                        border,
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)
summary(model1.4)

Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, 
   ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
   scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
   ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
   verbose, ", " lincomb = lincomb, selection = selection, control.compute 
   = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, 
   control.fixed = control.fixed, ", " control.mode = control.mode, 
   control.expert = control.expert, ", " control.hazard = control.hazard, 
   control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " 
   control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
   ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
   num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, 
   working.directory = working.directory, ", " silent = silent, inla.mode 
   = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = 
   .parent.frame)") 
Time used:
    Pre = 10.4, Running = 3.59, Post = 0.438, Total = 14.4 
Fixed effects:
                           mean    sd 0.025quant 0.5quant 0.975quant   mode kld
density_high             -0.001 0.019     -0.038   -0.001      0.037 -0.001   0
density_low              -0.033 0.014     -0.060   -0.033     -0.007 -0.033   0
sep1                      0.061 0.021      0.021    0.061      0.101  0.061   0
sep2                     -0.007 0.017     -0.041   -0.007      0.027 -0.007   0
sep4                     -0.030 0.017     -0.063   -0.030      0.002 -0.030   0
sep5                     -0.021 0.018     -0.056   -0.021      0.014 -0.021   0
border                    0.028 0.022     -0.015    0.028      0.072  0.029   0
sexFemale:age_group40-59  0.013 0.027     -0.041    0.013      0.067  0.013   0
sexMale:age_group40-59    0.243 0.023      0.199    0.243      0.287  0.243   0
sexFemale:age_group60-69  0.046 0.024      0.000    0.046      0.092  0.046   0
 [ reached getOption("max.print") -- omitted 5 rows ]

Random effects:
  Name    Model
    id_space BYM2 model

Model hyperparameters:
                         mean    sd 0.025quant 0.5quant 0.975quant  mode
Precision for id_space 64.272 8.855     48.954   63.534     83.767 61.89
Phi for id_space        0.527 0.102      0.331    0.526      0.728  0.52

Watanabe-Akaike information criterion (WAIC) ...: 93561.17
Effective number of parameters .................: 613.54

Marginal log-Likelihood:  -46519.56 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.4$summary.fixed)

                              mean       sd 0.025quant  0.5quant 0.975quant
density_high             0.9993970 1.019164  0.9628487 0.9993981  1.0373253
density_low              0.9674046 1.013621  0.9420980 0.9673945  0.9934492
sep1                     1.0627313 1.020778  1.0207214 1.0627293  1.1064822
sep2                     0.9927683 1.017501  0.9595510 0.9927680  1.0271377
sep4                     0.9700724 1.016812  0.9388166 0.9700878  1.0022794
sep5                     0.9791846 1.018250  0.9450946 0.9791696  1.0145921
border                   1.0288841 1.022512  0.9848690 1.0288992  1.0747769
sexFemale:age_group40-59 1.0128943 1.027798  0.9598695 1.0128928  1.0688572
sexMale:age_group40-59   1.2746180 1.022762  1.2195962 1.2746116  1.3321598
sexFemale:age_group60-69 1.0470351 1.023857  0.9997301 1.0470326  1.0965928
                              mode kld
density_high             0.9993998   1
density_low              0.9673744   1
sep1                     1.0627253   1
sep2                     0.9927673   1
sep4                     0.9701189   1
sep5                     0.9791396   1
border                   1.0289297   1
sexFemale:age_group40-59 1.0128897   1
sexMale:age_group40-59   1.2745988   1
sexFemale:age_group60-69 1.0470278   1
 [ reached 'max' / getOption("max.print") -- omitted 5 rows ]

model1.4$waic$waic - model1.3$waic$waic 

[1] -8.908195

drivers_plot(model1.4,data1)

We see that the model fit improves (lower WAIC) and that some of the spatial variability is not explained by the covariates (precision for id_space increases). We find two clear effects. First, low population density is associated with lower relative excess mortality. Second, the lowest quintile of socio-economic position (the 20% least affluent municipalities according to this metric) have higher relative excess mortality. In both these cases, the 95% credible intervals exclude 1, indicating a clear effect, but there are also other associations that may exist with border municipalities (in favour of a positive association) and with highest SEP quintiles Q4 and Q5 (in favour of a negative association).

All of these estimates apply similarly to all age and sex groups (no interaction), and are adjusted for the regional variability associated with COVID-19 waves in 2020.

Model 1.5: language regions

One idea is to look at the effect of language regions. The difficulty is the colinearity between language regions and the first COVID-19 wave of 2020, that primarily affected Ticino (Italian) and Southwestern Switzerland (French), mostly because of how the initial global spread of COVID-19 occurred (with large early epidemics in Italy and then France). These effects are much larger than any effect that could be attributed to cultural differences between language regions, so it is quite difficult to estimate the latter. We still attempt to do so by adding the language of each municipality (reference is German) to our model.

model1.5 = INLA::inla(munici_observed ~ - 1 + offset(E) +
                        sex:age_group +
                        f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                          hyper = hyper.bym2, constr=TRUE) +
                        density_high + density_low +
                        sep1+sep2+sep4+sep5 +
                        border +
                        lang_fr + lang_it,
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)
summary(model1.5)

Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, 
   ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
   scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
   ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
   verbose, ", " lincomb = lincomb, selection = selection, control.compute 
   = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, 
   control.fixed = control.fixed, ", " control.mode = control.mode, 
   control.expert = control.expert, ", " control.hazard = control.hazard, 
   control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " 
   control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
   ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
   num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, 
   working.directory = working.directory, ", " silent = silent, inla.mode 
   = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = 
   .parent.frame)") 
Time used:
    Pre = 10.2, Running = 3.86, Post = 0.442, Total = 14.5 
Fixed effects:
                           mean    sd 0.025quant 0.5quant 0.975quant   mode kld
density_high              0.004 0.019     -0.033    0.004      0.041  0.004   0
density_low              -0.030 0.014     -0.057   -0.030     -0.004 -0.030   0
sep1                      0.050 0.021      0.010    0.050      0.091  0.050   0
sep2                     -0.013 0.017     -0.047   -0.013      0.021 -0.013   0
sep4                     -0.028 0.017     -0.060   -0.028      0.005 -0.028   0
sep5                     -0.019 0.018     -0.054   -0.019      0.016 -0.019   0
border                    0.013 0.023     -0.031    0.013      0.058  0.013   0
lang_fr                   0.086 0.026      0.035    0.086      0.138  0.086   0
lang_it                   0.126 0.047      0.035    0.126      0.218  0.126   0
sexFemale:age_group40-59 -0.016 0.029     -0.072   -0.016      0.040 -0.016   0
 [ reached getOption("max.print") -- omitted 7 rows ]

Random effects:
  Name    Model
    id_space BYM2 model

Model hyperparameters:
                         mean    sd 0.025quant 0.5quant 0.975quant   mode
Precision for id_space 67.176 9.457     50.863   66.371     88.038 64.589
Phi for id_space        0.497 0.106      0.296    0.495      0.707  0.485

Watanabe-Akaike information criterion (WAIC) ...: 93551.93
Effective number of parameters .................: 608.22

Marginal log-Likelihood:  -46525.05 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.5$summary.fixed)

                              mean       sd 0.025quant  0.5quant 0.975quant
density_high             1.0037454 1.019052  0.9672432 1.0037476   1.041611
density_low              0.9701207 1.013597  0.9447862 0.9701108   0.996192
sep1                     1.0515420 1.020870  1.0097970 1.0515396   1.095027
sep2                     0.9867008 1.017509  0.9536698 0.9867006   1.020877
sep4                     0.9725319 1.016759  0.9412940 0.9725468   1.004720
sep5                     0.9810411 1.018143  0.9470841 0.9810248   1.016311
border                   1.0134256 1.022794  0.9696052 1.0134191   1.059265
lang_fr                  1.0902911 1.026380  1.0360236 1.0902752   1.147496
lang_it                  1.1344956 1.047944  1.0351661 1.1343880   1.244022
sexFemale:age_group40-59 0.9841035 1.028953  0.9305179 0.9841080   1.040747
                              mode kld
density_high             1.0037515   1
density_low              0.9700910   1
sep1                     1.0515348   1
sep2                     0.9867003   1
sep4                     0.9725768   1
sep5                     0.9809922   1
border                   1.0134053   1
lang_fr                  1.0902442   1
lang_it                  1.1341726   1
sexFemale:age_group40-59 0.9841172   1
 [ reached 'max' / getOption("max.print") -- omitted 7 rows ]

model1.5$waic$waic - model1.4$waic$waic 

[1] -9.235079

drivers_plot(model1.5,data1)

At first sight, we may thing that there is a large effect of language region on excess mortality, with 9% more deaths than expected in French-speaking municipalities and 13% more in Italian-speaking municipalities compared to German. However, as expected this association is likely counfounded by the regional variability associated with COVID-19 waves in 2020. Indeed, if we now look at the geographically-structured municipality effect, we see that the higher excess in South and Southwestern Switzerland is now more evenly distributed (captured by the language effect), while French-speaking regions that were not as impacted during the first wave (such as Neuchatel and Jura) now have a negative municipality effect to compensate. These nonsensical results highlight the difficulty to estimate the effect of language regions, which remains inconclusive.

Model 1.6: voting

We leave aside the question of language, and now focus on results from two referendums about COVID-19 control measures held in June and November 2020. The point here is not to look at causality one way or the other, as we look at overall excess for 2020, and the voting took place at two separated points. A preliminary analysis has reported a negative association between the proportion of “yes” vote at the November referendum at the cantonal level and 7-day incidence on December 7, 2021 (https://smw.ch/index.php/smw/announcement/view/50).

model1.6 = INLA::inla(munici_observed ~ - 1 + offset(E) +
                        sex:age_group +
                        f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                          hyper = hyper.bym2, constr=TRUE) +
                        density_high + density_low +
                        sep1+sep2+sep4+sep5 +
                        border +
                        vote_jun_q1 + vote_jun_q2 + vote_jun_q4 + vote_jun_q5,
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)
summary(model1.6)

Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, 
   ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
   scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
   ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
   verbose, ", " lincomb = lincomb, selection = selection, control.compute 
   = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, 
   control.fixed = control.fixed, ", " control.mode = control.mode, 
   control.expert = control.expert, ", " control.hazard = control.hazard, 
   control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " 
   control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
   ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
   num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, 
   working.directory = working.directory, ", " silent = silent, inla.mode 
   = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = 
   .parent.frame)") 
Time used:
    Pre = 10.4, Running = 4.09, Post = 0.456, Total = 14.9 
Fixed effects:
                           mean    sd 0.025quant 0.5quant 0.975quant   mode kld
density_high             -0.003 0.020     -0.042   -0.003      0.035 -0.003   0
density_low              -0.034 0.014     -0.062   -0.034     -0.006 -0.034   0
sep1                      0.063 0.021      0.023    0.063      0.104  0.063   0
sep2                     -0.005 0.018     -0.039   -0.005      0.030 -0.005   0
sep4                     -0.030 0.017     -0.063   -0.030      0.003 -0.030   0
sep5                     -0.022 0.019     -0.059   -0.022      0.015 -0.022   0
border                    0.028 0.022     -0.015    0.028      0.072  0.028   0
vote_jun_q1               0.014 0.021     -0.026    0.014      0.055  0.014   0
vote_jun_q2               0.036 0.018      0.000    0.036      0.072  0.036   0
vote_jun_q4               0.023 0.017     -0.010    0.023      0.056  0.023   0
 [ reached getOption("max.print") -- omitted 9 rows ]

Random effects:
  Name    Model
    id_space BYM2 model

Model hyperparameters:
                         mean    sd 0.025quant 0.5quant 0.975quant   mode
Precision for id_space 64.144 8.844     48.808   63.420     83.581 61.828
Phi for id_space        0.521 0.102      0.325    0.519      0.721  0.512

Watanabe-Akaike information criterion (WAIC) ...: 93564.96
Effective number of parameters .................: 619.50

Marginal log-Likelihood:  -46547.43 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.6$summary.fixed)

                  mean       sd 0.025quant  0.5quant 0.975quant      mode kld
density_high 0.9969608 1.019799  0.9593352 0.9969597  1.0360682 0.9969568   1
density_low  0.9666332 1.014381  0.9399654 0.9666224  0.9941206 0.9666007   1
sep1         1.0654565 1.020926  1.0230496 1.0654538  1.1096368 1.0654485   1
sep2         0.9954325 1.017674  0.9618061 0.9954313  1.0302419 0.9954289   1
sep4         0.9703819 1.017142  0.9385242 0.9703952  1.0032435 0.9704220   1
sep5         0.9778523 1.018983  0.9424809 0.9778361  1.0146464 0.9778036   1
border       1.0288405 1.022537  0.9847802 1.0288554  1.0747842 1.0288854   1
vote_jun_q1  1.0145749 1.020944  0.9741284 1.0145840  1.0566472 1.0146022   1
vote_jun_q2  1.0363657 1.018452  0.9998833 1.0363561  1.0742356 1.0363367   1
vote_jun_q4  1.0232151 1.016912  0.9901118 1.0232110  1.0574488 1.0232029   1
 [ reached 'max' / getOption("max.print") -- omitted 9 rows ]

model1.6$waic$waic - model1.4$waic$waic 

[1] 3.793261

drivers_plot(model1.6,data1)

The results are not conclusive. Compared to the 3rd quintile of “yes” vote to maintain the COVID law at the June referendum, taken as a reference, we observe slightly increased relative excess mortality in all other quintiles, always compatible with no effect (relative excess = 1). This only suggests a lower excess mortality in the third quantile, corresponding to municipalities with between 53 and 65% of “yes”, which is nonsensical. This can be shown by changing the reference group for instance to Q1.

model1.6b = INLA::inla(munici_observed ~ - 1 + offset(E) +
                        sex:age_group +
                        f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                          hyper = hyper.bym2, constr=TRUE) +
                        density_high + density_low +
                        sep1+sep2+sep4+sep5 +
                        border +
                        vote_jun_q2 + vote_jun_q3 + vote_jun_q4 + vote_jun_q5,
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)

drivers_plot2(model1.6b,data1)

Taking Q1 as a reference highlights the absence of clear association between results of the June referendum and 2020 excess mortality in this analysis.

model1.6c = INLA::inla(munici_observed ~ - 1 + offset(E) +
                        sex:age_group +
                        f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                          hyper = hyper.bym2, constr=TRUE) +
                        density_high + density_low +
                        sep1+sep2+sep4+sep5 +
                        border +
                        vote_nov_q1 + vote_nov_q2 + vote_nov_q4 + vote_nov_q5,
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)

drivers_plot(model1.6c,data1)

We also don’t find any clear association when looking at the November referendum.